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Abstract 

A two-component model is developed consisting of a discrete loop of cardiac cells that circu- 
lates action potentials as well as a pacing mechanism. Physiological properties of cells such 
as restitutions of refractoriness and of conduction velocity are given via experimentally mea- 
sured functions. The dynamics of circulating pulses and their interactions with the pacer are 
regulated by two threshold relations. Patterns of spontaneous initiations and terminations of 
reentry (SITR) generated by this system are studied through numerical simulations and an- 
alytical observations. These patterns can be regular or irregular; causes of irregularities are 
identified as the threshold bistability of reentrant circulation (T-bistability) and in some cases, 
also phase-resetting interactions with the pacer. 
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1 Introduction. 

Ventricular arrhythmia is the leading cause of cardiac arrest and sudden death. Clinical observa- 
tions and implantable cardioverter defibrillators have accumulated a substantial amount of data 
on the occurrences of ventricular arrhythmia in patients [1, 36, 37, 38, 52, 59-61]. Temporal pat- 
terns of initiations and terminations of ventricular arrhythmia tend to exhibit substantial variations 
across different time scales, and their occurrences do not correlate decisively with medication, exer- 
tion, stress, lifestyles and similar factors. Arrhythmia events are not random; they show circadian 
patterns [36, 60] and also tend to occur in clusters. However, the detection times between consec- 
utive events and clusters are spread out over time [37, 52, 59], making it difficult to understand 
their causes and make predictions about their occurrences. Unlike the circadian patterns, these 
clusterings or their patterns of occurrences are not affected by the long-term administration of 
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antiarrhythmic drugs [59]. In spite of the abundant data in existence, the basis for non-circadian 
patterns is not well- understood. 

Many factors, ranging from internal cardiac mechanisms to external chance events play signif- 
icant roles in shaping the electrocardiogram (ECG) recordings and the implantable defibrillator 
data. Separating all of the possible contributions is a formidable task, but understanding the influ- 
ences of various factors and the extent to which each plays a role may have important consequences 
for the diagnosis and treatment of tachyarrhythmias. A number of prior studies using special prepa- 
rations of animal cell cultures establish that complex, spontaneously generated patterns may occur 
without some of the features peculiar to the heart (e.g. the 3-dimensional geometry or biological 
features such as valves, different tissue layers and types, etc) [3, 25, 35, 44]. These studies also 
indicate (both theoretically and experimentally) that causes external to the heart, whether random 
or deterministic, are not always necessary for the occurrence of irregular behavior. A different 
class of studies that involve interactive self-oscillatory sources also establish the capability of basic 
cardiac mechanisms to generate complex rhythms without considering the whole heart [23, 28, 43]. 

From these and similar studies it may be inferred that changes in the basic internal cardiac 
mechanisms can play a significant role in generating and sustaining arrhythmias. To clarify the 
role of these basic internal mechanisms at a fundamental level, in this paper we focus on the unidi- 
rectional reentrant circulation of action potentials in a loop of ventricular tissue. Because of their 
pervasiveness, relative simplicity and importance in clinical arrhythmogenesis reentrant loops have 
received a great deal of attention [4, 8, 11, 14-16, 20, 26, 30, 41, 44, 48-50, 55]. Studying spon- 
taneous initiation and termination of reentry (SITR) patterns in the loop context is considerably 
simpler than in the whole heart and can offer potentially useful insights into the clinically observed 
patterns of arrhythmia. In the construction of a workable model, a certain level of abstraction 
bridges over physiological complexities that do not play a central role in the long-term evolution of 
temporal patterns. 

In this paper we consider a model that combines the traditional discrete loop with a pacer to 
form an interactive system. Fundamental thresholds are dynamically integrated so they can either 
trigger reentry or inhibit it. Although a relatively simple mechanism, this composite system has 
two important features: It is an essential component of a well-known cardiac anomaly; and it is 
capable of generating complex initiation and termination patterns. 

In a series of case-studies we show that the loop-pacer mechanism may be responsible for irreg- 
ularities in onset and termination of reentry beyond such pervasive factors as tissue heterogeneity, 
the geometry of the heart and other physiological and non-physiological factors. Even without 
such complexities, the SITR patterns generated by the loop-pacer system can be complex; i.e. they 
contain many regular features yet their evolution over time is difficult to predict. We trace this 
complexity to the interaction of a threshold with bistability of the reentrant circulation and in 
some cases, also to the inherent discontinuity of phase resetting interactions between the loop and 
the pacer. Bistability is an inherent physiological characteristic that may be attributed to non- 
uniformities in the restitution of action potential duration (APD); this feature emerges prominently 
when conduction velocity (CV) is properly matched with APD. When this happens (e.g. when the 
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Figure 1: Schematic diagram of a loop joined to the rest of the heart by a single conducting pathway. 

length of the loop is within a certain range) the SITR patterns become considerably more complex 
and unpredictable. 

2 The threshold model. 

In this section we describe the main features of the loop-pacer model, leaving a few additional 
details to Section 4. The pacer is a self-oscillatory mechanism that is tied to the reentrant circuit 
(the loop in this paper) via a time lag or delay parameter and a special threshold. 

2.1 The loop. 

Consider a loop of cardiac tissue consisting of cells that conduct action potentials and assume that 
a single conducting pathway connects the loop to the rest of the heart. Let the physical length of 
the loop be denoted by L measured in centimeters. Divide the loop into m sets or aggregates of 
cells, each of which may be called a cell aggregate or for brevity, just a "cell." The cell aggregate 
that is connected to the only pathway out of the loop is also the gateway through which pulses 
enter the loop or exit it; we label it Cell 1; see Figure 1. This cell is adjacent not only to Cell 2, 
but also to Cell m at its other end. 

For each integer i between 1 and m, label the length of the i-th Cell ALj. These cell aggregates' 
lengths are not necessarily equal but of course, they add up to L. Since each cell aggregate must 
have at least one cardiac cell in it, ALj has a natural lower bound, namely the nominal length of 
a single cardiac cell (about 0.01 cm or 100 microns). The closer all ALj are to this lower bound, 
the greater is the number m of cell aggregates in a loop of fixed length L; for a 12 cm loop with 
cells of length 0.01 cm, the number of cells is 1200. While for discrete modeling it is not necessary 
to pick m that large, it needs to be large enough (hence each cell aggregate small enough) that the 
conduction velocity from one end of an aggregate to the other can be taken to be approximately 
constant. This assumption is essential if the restitution of conduction time (CT) is defined in terms 
of the restitution of CV, as we do later on. 

We assume that all the cardiac cells within a given aggregate are identical. If all cells in 
the entire loop have identical physiological characteristics then we may also let cell aggregates all 
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have equal lengths. In such a case, the loop is said to be "homogeneous." Otherwise the loop is 
"heterogeneous ." 

To initiate unidirectional circulation, it is necessary that a unidirectional block (UB) exist 
somewhere in the loop; for simplicity, we place it in Cell 1, the gateway to the loop; see Figure 
1. Thus Cell 1 is able to activate Cell 2 but not Cell m. For our simulations, we use a simplified 
version of the UB time window that is defined by a threshold value for the diastolic interval (DI). 
Thus conduction is blocked in any cell whose DI is not greater than the threshold DI* ; see Section 
4 below for more details. Such a minimum value is used in [26] as a termination mechanism for 
explaining the experimental results of [20]; also see [19]. 

2.2 Restitution functions. 

We use the term "restitution" generally for relations that give a particular quantity, e.g. action 
potential duration or conduction time, as a function of the diastolic interval (DI). The DI (usually 
measured in milliseconds, ms) is the rest or recovery period for the cell which essentially starts with 
the end of the refractory period and ends when the cell is activated again. 

2.2.1 APD restitution. 

In its most basic form, the action potential duration or APD is the length of time (usually measured 
in milliseconds, ms) that a cell is active after excitation. For our purposes, we may think of APD as 
a cell's effective refractory period (ERP) during which no excitations, even strong ones, can elicit 
new action potentials. The restitution of APD, which can be experimentally measured or derived 
from ionic models, is the most extensively studied of restitution relations [2, 4-8, 10, 11, 21, 24, 26, 
29, 32, 33, 42, 45, 48, 50, 53-57]. Any APD restitution function, whether experimentally measured 
or analytically derived can be used in this model. For the numerical simulations in this paper, we 
use an APD restitution function that is fitted to the experimental data by Koller, et al [33] where 
APD values were recorded in two types of patients, those with and those without structural heart 
disease (SHD). The following is a possible fit to their averaged data displayed in Figure 1 in [33] 
for SHD patients: 

A(DI) = ai - a 2 e-^ DI - a^ DI -^ - a 4 e^(^- 2 ) 2 + ^(DI- t 3 ) 

(DI - T 3 y + a 6 

with parameter values 
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A graph of this function is shown in Figure 2. Expressing the APD as a function of DI in the 
above manner represents only a first order of approximation (see Section 4 below). 
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Figure 2: The APD restitution curve given by fTJ). The dashed line has slope 1 and is included for reference. 
2.2.2 CT restitution. 

Each cell conducts an action potential through it in a finite amount of time. This time interval, 
usually measured in milliseconds, is referred to as the conduction time (CT) [4, 26, 50]. Experi- 
mentally measured CT restitution functions are not as readily available for human hearts as the 
APD restitution functions; however, they may be readily derived from CV restitution functions 
through the relation 

C < D " = vm < 2 > 

where C and V are, respectively, the CT and the CV restitution functions. This derivation is valid 
if as mentioned earlier, V is essentially constant over the length of each cell (or cell aggregate) and 
thus no spatial dependence is required in V. 

For discussions of CV restitution functions see [2, 7, 8, 11, 12, 14-16, 20, 22, 53, 54, 56]. For 
our simulations we use the following expression 

C(DI) = ^ [1 + de-^ 1 ] , DI > (3) 

with parameter values: 



AL 


c 


d 




0.1 cm 


0.07 cm/ms 


1 


0.02 ms" 1 



These numbers are obtained using ([2j) and a reasonable fit to the human-rescaled, guinea pig 
data for CV from [22, 54]. The value of AL corresponds to 10 nominal cells. The number c in ([3j)is 
the maximum CV; in more common (though less uniform) units, the fitted number in the above 
table would be 0.7 m/s. By changing the parameters in (J3j) we may achieve conduction slow-downs 
in different ways, as shown in Figure 3. 
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Figure 3: CT restitution curves given by (O for several choices of parameters. For each curve, AL = 0.1 
and c = 0.07. (a) d = 1 and w = 0.02; (b) d = 1.5 and uj = 0.02; and (c) d = 1 and lo = 0.01. 

2.3 The two rhythms. 

In the absence of reentry, a pacer (the sinus node or an ectopic focus) drives the ventricular 
contraction cycle. If reentry is initiated and interrupted then the pacer and the loop fall out of 
phase and it is necessary to carefully model their interaction. For this purpose we define two 
rhythms or clocks: The pacer's rhythm which has a prescribed beat pattern (fixed or variable) and 
the system's rhythm which is reset when reentry is initiated or terminated. These two rhythms are 
related via a pacer threshold (seen the next section) and a delay term that we now define. 

The time of occurrence for the n-th beat (following a fixed reference beat) can be tracked and 
calculated using the following quantity, namely, the sum of usually variable cycle lengths (after a 
DI adjustment): 



Let the time of occurrence of the k-th pacer beat be denoted by (starting from the same 
reference point in time as that for p n ). The difference 



is the duration of each such beat. In our simulations B(k) is defined arbitrarily in order to facilitate 
our study of the SITR patterns. If the pacer is the sinus node then B(k) is influenced by a variety 
of deterministic and stochastic factors, such as circadian, weekly and seasonal variations, lifestyles, 
drugs, the autonomic system and so on. Thus a single, all purpose formula for B(k) does not exist. 
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Pn = CLj + (DI lfi - DI 1>n ), CL, = A(Z)/ lti _i) + DI hj . 



Ph ~ 0k-i = B{k) 



G 



The two rhythms p n and (3k fall out of phase if reentry occurs as noted earlier. To model their 
interaction, we define the first pacer pulse that reaches Cell 1 in beat n as the least integer k n such 
that 

Pk n >pn + S n . (4) 

The variable delay term 5 n here represents a time lag that results from the blocking of the pacer 
by the retrograde reentrant wave propagating toward the pacer. A number of things may affect <5 n ; 
these include the different types of tissue in which pulses propagate (retrograde and antegrade), 
the electrotonic spread of current in conducting tissue between the pacer and the loop, and the 
effect of reentrant waves on the pacer's intrinsic beat rate. For simplicity we do not consider the 
electrotonic currents and tissue heterogeneity. We assume also that the pacer itself is protected in 
the sense that it's intrinsic beat rate is not affected by the reentrant activity. 

Therefore, we use a simple definition for 5 n as follows: 

5, if pulse n — 1 was reentrant ^ , 
DI*, otherwise 

If the loop is not in reentry, then 6 n may take the minimum value DI* since Cell 1 must have 
at least that much rest time after its ERP before it can be reactivated. High frequency reentrant 
waves typically disassociate a slow pacer such as the sinus node. This situation is modeled here 
by choosing a sufficiently large value of the delay parameter 5 so as to inhibit the pacer from 
interfering with the reentrant circulation. Small values for 5, which promote pacer interference, 
may be feasible in certain circumstances, e.g. when there are pathways from the pacer's site to the 
site of the loop that are protected from the fast reentrant waves by pathologies. 

To determine the index value k n we may generate the sequence independently of p n and 
for each n mark the values k n that satisfy (|4]). This process is carried out internally during each 
simulation run; an analytical relationship is generally not easy to find. However, if B(k) = Bq is 
constant then such a relationship is easy to find; in this case we have (3k n = k n Bo with 



kn 



Pn + &n 

B n 



i.e., the least integer that is greater than or equal to the ratio (p n + S n )/ Bq. This simple relationship 
makes explicit the discontinuity that results from phase resetting in our model. 



2.4 Modes and thresholds. 

We distinguish between two primary dynamic modes for the loop-pacer system: The reentry mode 
where an action potential circulates in the loop by itself and the paced mode where the action 
potential in the loop comes from the pacer. The system's mode in each beat changes on the 
basis of threshold relations to be defined below. The modes are defined by systems of partial 
difference equations that track the circulating pulse in both space and time. Introductory material 
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on ordinary and partial difference equations can be found in [9, 17, 31, 40, 51]. Technically, the 
dynamical system consisting of the loop, the pacer and the associated thresholds is a polymodal 
structure in the sense of [51]. 



2.4.1 Dynamic activation duration. 

Let DIi n be the DI of Cell i in beat n. Define the firing indicator function as 

0i,n = 1) and 

1, if Dli n _i > DI* and 0;_i >n = 1 



i = Til 

0, if DJ iin _! < DI* or &_ 1>n = ' Z ' ' ' ' ' m ' 

This quantity specifies conditions under which the i-th cell fires an action potential in beat n; 
these conditions simply require that (a) the preceding Cell i — 1 fired and (b) the DI of Cell i in 
the preceding beat n — 1 exceeded the minimum value DI* as is required for propagation. 

We assume that Cell 1 fires in every beat whether by a reentrant pulse or by a pulse from the 
pacer. 

Let APDi n denote the action potential duration of Cell % in beat n so that if the i-ih cell does 
not fire in beat n then APDi^ n = (otherwise, APDi^ n = Aj(Z)Jj jn _x))- The dynamic activation 
duration is defined as 

AD i>n = 4>i tn Ai(DI i)n -i), i = 1, 2, . . . , m 

where Ai is the APD restitution function for Cell i. Note that AD\^ n = APD\^ n = Ai(DJi )n _i) for 
Cell 1 in every beat. In the following sections (except in Case-study 4) we assume that the loop is 
APD homogeneous so that A{ = A for alii = 1, 2, . . . , m where A is given by (pQ). 



2.4.2 Paced mode. 

In this mode the pacer determines the DI values DI{ n rather than a reentrant pulse in the primary 
loop; the circulation is assumed unidirectional owing to an active UB. We use the following system 
of m equations: 



DIi >n = P kn ~ Pn (5a) 
DJ ilB = £>!<_!,„ + AA-i,n - AD <iB + CT^^ - CT^ 1)n (5b) 

+ (1 - 4>i,n)DIi,n-l ~ (1 - <t>i-l,n)DIi-l,n-l 

with i = 2, . . . , m. 

Here C{ is the CT restitution function for Cell i. The above system of equations determines the 
DI values using the first pacer beat after time instance p n (beat n in progress). Figure 4 illustrates 
equations (|5bp and (|6bp . The first five terms in (|5b|) simply express in mathematical terms what 
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Figure 4: Schematic diagram of a circulating action potential in a ring of m cells; A is APD, D is 
DI and CT is the conduction time. The sizes of A, DI and CT are selected arbitrarily for clarity of 
illustration. 

we might state less precisely using the English language. They are the same for tissue fibers as 
they are for loops; see, e.g. [68]. The last two terms in the definition of DIi :U in (|5b[) allow for 
proper adjustment of the DI value (without double counting) when Cell i does not fire in beat n 
(see comments on the reentry mode below). In typical loop models in prior literature all cells fire 
so these two terms would drop out. 

2.4.3 Reentry mode. 

In this mode the reentrant loop drives the ventricles ectopically. Thus Dl\^ n is the difference 
between the APD of Cell 1 and the total conduction time for unidirectional circulation once around 
the ring, i.e. 

m 

Dh, n = J2 CT^ n - A\(DI\ tn -\) (6a) 

3=1 

DI hn = Dh_ 1)n + AA_i >n " AD ijW + CTi_ 1>n+1 - CTi- ltn (6b) 

+ (1 - 0i,n)DIi t n-i - (1 - 0i_i )n )jDii_i j?l _i 

with % = 2, . . . , m. 

It is not hard to show that the system of m equations (J6j) is equivalent to the coupled map 
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lattice equations in [51] when = 1 for all i and n. 

Equations ([6]) are not valid if some of the DI values fall below the threshold DI* or if a stray 
pacer pulse gets into the loop (i.e. if the right hand side of (foa|) exceeds the right hand side of 
(|5aj) ). The following threshold relations ensure a consistent scenario. 

TC: The circulation threshold (or head-tail or conduction block threshold) is based on 
the following set of inequalities: 

in 

Y,CT jin -A l {DI l ^ 1 )>DI*, 

3=1 

DIi-i,n + ADi-. x ,n - AD ijn + CTj_i >n+ i - CT^n > DI* for i = 2, . . . , m. 

They ensure that the DI values DIi^ n in Equations © exceed DI*, as is necessary for 
the generation of action potentials. If any of the above inequalities fails in a Cell i then 
a wavefront reaches its own tail in that particular cell and circulation stops until the 
arrival of the next pacer pulse. In this case we say that TC fails in Cell i. 



TP: The pacer threshold (or phase resetting threshold) is defined by the following in- 
equality: 

m 

CT hn - A x {DIx, n -i) < f3 kn - Pn . 

3=1 

This inequality ensures that the reentrant pulse (rather than a stray one from the pacer) 
reactivates Cell 1. 

Mode switching criteria. Equations (GJ) are used if TP holds and also TC holds in Cell 1. 
If for some Cell i > 2 TC fails, then (I6bj) can still be used to find the DI values for beat n, but we 
switch to equations ^ for beat n + 1. If TC fails for i = 1, or if TP fails then we use ([5]j for beat 
n. 

2.4.4 Bidirectional circulation (BDC) mode. 

In the absence of unidirectional block the pacer sends two pulses in the loop with Cell 1 activating 
both Cell 2 and Cell m. These pulses propagate in opposite directions to eventually annihilate each 
other in some cell in the loop. Transitions into and out of this bidirectional circulation mode are 
governed by TC, TP and a UB threshold (TUB); see Section 4 below. We assume for simplicity 
that TUB never fails so that the BDC mode need not be discussed in this paper. 
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3 SITR patterns. 



In this section we present the results of numerical simulations of our threshold model along with 
some analytical observations. Due to the large number of different types of behavior that can 
occur by changing one or more of the many variables involved, providing a comprehensive list of 
essentially different patterns is not feasible here. We therefore present our main results in the form 
of a few model case-studies and for additional results we refer the interested reader to our web site: 
www.people.vcu.edu/~hsedagha/SITR. 

Each case-study below consists of one or more "runs," i.e. sets of iterations of the propagation 
equations. Each run is based on a fixed set of parameter values and we call each iteration in a 
given run a "beat". For each run we need to specify m initial DI values Dlip. This specification of 
DI values is a technical representation of the occurrence of a premature stimulation in our model. 

Computational accuracy is not a major issue in the qualitative studies in this paper. However, 
calculating the values of the various nonlinear functions (e.g. the exponentials in the APD and CV 
restitution functions) generate round-off errors that can accumulate in certain cases (e.g. runs in 
Case-study 4). Therefore, results obtained by different programs, codes or levels of precision (single 
or double) may be different from some of those that we have presented, although the qualitative 
features are usually retained. 

3.1 Mode sequences. 

The following conventions facilitate the presentation of results in this section: 

For a loop consisting of m cells if a reentrant pulse is blocked at Cell j in a particular 
beat, then we define the reentry value of that beat as the ratio (j — l)/m; this ratio will 
be abbreviated as [j]. 

If a reentrant pulse completes a turn around the loop then the reentry value is 1. 

If a reentrant pulse is blocked in Cell j after t complete turns around the loop then the 
reentry value is t[j]. 

If a beat occurs in the paced mode (i.e. one of the thresholds TC or TP fails) then we 
assign it the value —1. Thus t consecutive paced beats appear as —t. 

If the pattern is locked in the reentry mode, we use the symbol oo; if it is locked in the 
paced mode, we use — oo. 

Adding up consecutive reentry mode beats and consecutive paced mode beats gives a sequence 
of numbers consisting of integers whose sign tell us the mode of system in various sets of beats. 
Examples of mode sequences and more details on them are found in the case-studies below. 
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3.2 Bistability and thresholds. 

We call the loop-pacer system bistable if there are two or more coexisting stable DI configurations 
or state vectors, each of which can be reached from a particular region of the m-dimensional state 
space. If there are more than two distinct, coexisting stable states then the system is multistable. 

Bistability may be attributed to the non-concavity of the APD restitution function [50]. The 
APD function ([1]) whose graph has bumps and twists is also non-concave. Proper CT restitution 
parameters are required to realize bistability; i.e. the CT and APD restitution parameters must 
be properly matched. Because the length L of the loop affects the conduction time through it, 
bistability may emerge or fade as L is changed with other CT parameters fixed [50]. 

A bistable regime can affect SITR patterns by causing the violation of the circulation threshold 
TC without any changes in the APD or CT parameters, or any changes in the pacing rate. This 
situation occurs if the bistable regime satisfies the following threshold bistability condition: 

T -bistability: There are two distinct stable states in the reentry mode: In one state the 
DI values cross DI* and cause the failure of the circulation threshold TC at some cell 
within the loop, but in the other state the DI values are always greater than DI* . 

An example of T-bistability is shown in Figure 7 below; other examples of T-bistability occur in 
the next section (Case-study 1). The m-dimensional state space of a T-bistable system is partitioned 
into two regions or basins of attraction. As a SITR pattern evolves with the motion of the DI state 
vector, mode changes occur if the state vector crosses over into the basin of attraction of a different 
stable regime. 

3.3 Case-study 1: T-bistable reentry. 

In this section we establish that T-bistability lends a measure of unpredictability to SITR patterns 
that is not due to phase resetting disruptions by the pacer (i.e. TP failures). For instance, com- 
paring the mode sequences in Runs 1.2, 1.4 and 1.6 below we see that the initial few bursts of fast, 
reentrant beats provide no obvious clues about the eventual modes that the patterns lock into. 

Consider a homogeneous loop with the APD and CT restitutions given by (pQ) and Q. We also 
set the length L and other parameters as: 



L 


m 


DI* 


5 


B 


12.5 cm 


125 


15.3 ms 


120 ms 


320 ms 



where Bq is the fixed cycle length of the pacer; it is chosen small in this case-study to facilitate re- 
initiations of reentry. The value of 5 is set high enough to inhibit pacer interference (TP never fails) . 
The following table summarizes the results of simulations; it is followed by a series of observations 
and elaborations. 
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Run No. 


DJ i>0 , i = 1, ... 125 


SITR pattern mode sequence 


1.1 


200 ms 


{-2,2,-2,2,...} 


1.2 


100 ms 


{-2,4,-2,32[73],-2,oo} 


1.3 


80 ms 


{-l,oo} 


1.4 


70 ms 


{4, -2, 14, -2,4,-2,22(49], -00} 


1.5 


60 ms 


{00} 


1.6 


50 ms 


{ll,-2,4,-2,27[57],-2,oo} 



(i) In Run 1.1 the SITR pattern immediately locks into a regular form of 2 paced beats followed 
by 2 reentrant beats. In Runs 1.2, 1.3, 1.5 and 1.6 the SITR pattern locks into the reentry mode. 
In Run 1.4 the paced mode is locked into. Note that the only thing that is changing from one run 
to the next is the initial DI values. 

(ii) In Runs 1.2, 1.3, 1.5 and 1.6 the DI values in the sustained or locked reentry modes converge 
to a fixed number and the eventual values of the fundamental parameters are 



DI 


APD 


cycle length CL = APD+DI 


59.4 ms 


173.6 ms 


233 ms 



The fixed, limiting DI value 59.4 represents a stable convergent state. This equilibrium DI is 
locally stable (attracting) because using the restitution functions ([1]) and (|3|) we compute [4, 26] 

A' (59.4) + C"(59.4) = 0.90 < 1. 

(iii) In Runs 1.2, 1.4 and 1.6 all the spontaneously terminated reentry bursts have oscillating 
(quasiperiodic) DI. Thus APD and CL also oscillate similarly. The oscillatory state is stable 
and eventually crosses DI* = 15.3 ms. Hence there is T-bistability (the oscillatory one and the 
convergent one) which is responsible for mode changes in each run. 

For example, in Run 1.2 the 2nd initiation of reentry (the 32-beat burst) starts with the DI 
state-vector in the basin of attraction of an oscillatory state which expands enough for the DI 
values to cross DI* (TC fails); see Figure 5. This terminates the second reentry burst at Cell 73 
in the 41st beat. This is the same type of termination mechanism as that in [19, 20, 26]; this type 
of oscillations in DI are also seen in [56] where they are related to the T wave alternans in ECG 
recordings; also see [35] in this regard. With the 3rd initiation of reentry, the DI state-vector falls 
in the basin of attraction of the convergent state which then effectively "traps" the DI trajectory 
and locks the SITR pattern into the reentry mode (here then, simplicity of behavior is not a good 
thing). 

(iv) Examination of the simulation data for Run 1.4 shows a recurrent failing of TC in Cell 1 
after the last burst of 23 reentrant beats. For all sufficiently large values of n, 
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Figure 5: Results of Run 1.2, showing variation of DI in Beats 10 through 50. The dashed line corresponds 
to DI* = 15.3. Amplification of oscillations in DI leads to termination of reentry in Beat 41, Cell 73 because 
-D^73,4i < DI*. With Beat 44 reentry is re-initiated and the system locks into the stable convergent state 
(after some transient oscillations). 

The number 96.6 is the equilibrium DI for the paced mode with pacing period Bq = 320. The 
paced-mode equilibrium happens to be stable because in the absence of reentry the paced mode is 
governed by the one-dimensional difference equation 

DI hn = F(Z)Ji, n _i) = B - A(D/i, n _i) 

with F'(96.6) = -A' (96.6) = -0.901; i.e. |F'(96.6)| < 1. Why the transition to the stable paced 
equilibrium occurs in this run and not others is unclear. Figure 6 shows the changes in DI values 
in this run. 

(iv) In Run 1.5, the given initial DI values put the state vector in the basin of attraction of the 
convergent state of the system. This is evidently not the case in Run 1.6 with even shorter DI. In 
Run 1.6, as in Run 1.4, the initial mode is oscillating-DI reentry. 

3.4 Case-study 2: Non-T-bistable reentry. 

We now consider some longer loops than in Case-study 1 for comparison. In these loops bistability 
is present but does not satisfy the T-bistability condition. A smaller variety of different SITR 
patterns are obtained (with all other system parameters having the same values as in Case-study 
1). We consider two different lengths: 

L = 13 cm (m = 130) where stable oscillating states exist all of which fail TC, but the 
convergent mode is absent since the DI fixed point of about 62.5 ms is unstable [4, 26] 

A'(62.5) + C"(62.5) = 1.013 > 1. 
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Figure 6: Results of Run 1.4, showing variation of DI in Beats 1 through 60. The dashed line corresponds 
to DI* = 15.3. The four interruptions of reentry due to conduction block (TC fails) are seen as brief jumps 
to high DI. The last interruption terminates reentry and locks the system in the paced mode. 

L = 14 cm (m = 140) where all of the stable oscillating states (there is at least one) 
satisfy TC. The convergent mode is again absent for the same reason as above: 

A'(68) + C'(68) = 1.253 > 1. 

To save space let us abbreviate repetitions or locked mode patterns with bars: 

k, -j = k, -j, k, -j, ... (7) 

The following table summarizes the results of simulations; it is followed by a series of observa- 
tions and elaborations. 



Run No. 


DIi, 


L = 13 cm 


L = 14 cm 


2.1 


200 ms 


{-2,2} 


{-2,2,-2,oo} 


2.2 


100 ms 


|-2,4,-2,12[60j,-2,2} 


{-2,oo} 


2.3 


80 ms 


{2,-2} 


{-l,l,-2,2,-2,oo} 


2.4 


70 ms 


{l2,-2,10[79],-2,17[79],-2,5[79]} 


{12,-2,2,-2,oo} 


2.5 


60 ms 


{56,-2,2} 


{-l,10,-2,9,-2,oo} 


2.6 


50 ms 


{3, -2,4,-2,12160], -2, 2) 


{-00} 



(i) The SITR patterns for the 14 cm loop are expected to lock into the reentry mode since the 
stable states do not fail TC to cause terminations. The exception is Run 2.6 where possibly owing 
to the greater length of the loop, a pulse from the pacer has time to activate Cell 1 and block the 



15 



premature simulation from reentering the loop (threshold TP fails). It is not clear why reentry 
does not re-initiate as in Run 2.5. 

(ii) In the case of the 13 cm loop the repeated terminations can be attributed to the fact that 
the stable states all fail TC; thus the patterns cannot lock into the reentry mode. 

(iii) For both of these loops, the eventual form of the SITR pattern is more predictable than 
the T-bistable loop of Case-study 1. On the other hand, the transient bursts of fast beats still seem 
difficult to predict. 

3.5 Case-study 3: Special patterns. 

Run 3.1 (more T-bistability). This run presents a complex SITR pattern that locks into a long 
cycle of several initiations and terminations. Make the following changes to parameters in Run 1.2: 



d 


B 


1.3 


315 ms 



Note that the increase in d from 1 to 1.3 in this run causes a slight elevation of the CT restitution 
curve (i.e., a conduction slow-down) at low DI values; see Figure 3. 

(i) In this run, after a transient period of spontaneous initiations and terminations, a long 62- 
beat pattern emerges that is repeated; i.e. the SITR pattern locks into a 62-beat cycle containing 
several bursts of fast reentrant beats. The mode sequence for this run is as follows, with the pattern 
between starred numbers repeating: 

{-2, 8[81], -1, [35], -2, 4, -2, 7[24], -1, [18], -2, 8[18],* -1, [12], -2, 7[12], -1, [12], -2, 
10[49], -2, 5[49], -2, 4, -2, 7[22], -1, [19], -2, 7[19],* -1, [12], -2, 7[12], -1, [12], -2, . . .} 

(ii) The DI equilibrium of approximately 64 ms is unstable in this run since 

A' (64) + C"(64) = 1.073 > 1. 

This instability is caused by the greater value of d and results in the absence of the convergent 
state. Nevertheless, there is T-bistability where one of the stable oscillating states does not fail 
TC but another oscillating state does. In the first oscillatory state, DI values come very close to 
DI* ; see Figure 7. This proximity increases the sensitivity to transient effects of mode changes and 
makes locking into the reentry mode an improbable event. 

Run 3.2 (pacer action). In Case-studies 1 and 2 we see that the convergent reentry pattern 
is "sticky" i.e. once this pattern is attained, reentry does not self-terminate by the circulation 
mechanism in the loop. However, pacer interaction with the loop can lead to termination in some 
cases by changing the convergent pattern to a terminating oscillatory one in a T-bistable case. 

All parameters are as in Run 1.5 (hence there is T-bistability) except that now: 
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Figure 7: Two distinct stable states illustrating T-bistability in Run 3.1. (a) If -DX^o = 67 for all i, the 
graph crosses the dashed horizontal line DI = 15.3; for the sake of illustrating the quasiperiodic nature of 
oscillations in reentry mode, the graph ignores TC since reentry would terminate with the first such crossing, 
(b) Dlifi = 66 for all i. Reentry is sustained with quasi-periodic variation of DI. 



5 


B 


45 ms 


800 ms 



We have set the pacing period at the nominal sinus length to eliminate the effects of fast pacing. 
The smaller value of 5 promotes interference by the pacer. 

(i) The mode sequence for this run, which may be compared with Run 1.5, is: 

{40,-l,15,-oo}. 

(ii) Reentry is initiated with Beat 1 by a premature stimulation, and the DI values begin to 
approach the limiting value 59.4 ms as in Run 1.5. However, TP fails in Beat 41; a stray pulse 
from the pacer enters the loop and changes the reentry DI pattern from the convergent type to an 
oscillating one; see Figure 8. The quantities PL n , PDI\ yn in Figure 8 are the right hand sides of 
equations ([5a]) and (f6a|) . respectively. 

The oscillating DI pattern crosses the DI* threshold (TC fails) in Beat 55 and reentry terminates 
in this beat. 



3.6 Case-study 4: Slow pacing, heterogeneity and PVC's. 

When the pacing period is long (say, within the nominal sinus range of 800±350 ms) and the loop 
is homogeneous, re-initiations of reentry do not occur with the above CT and APD parameters in 
the absence of premature stimulations or ectopic sources. The low CT (or high CV) in the loop 
blocks reentry when pacing is slow. Nevertheless, the rather common occurrences of ventricular 
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Figure 8: Results of Run 3.2. (a) Variations in PL n — PDl\ n (see text). The negative value in Beat 41 
(indicated by the arrow) corresponds to termination of reentry due to failure of TP. (b) DI values in Beats 
30 through 60; note the phase resetting that occurs in Beat 41 and changes the reentrant circulation pattern 
from convergent to oscillatory in this T-bistable case. 

tachycardia at low pacing rates in patients [58] motivate us to study this important special case in 
the context of a heterogeneous loop with a variable-rate pacer. 

Suppose that the loop consists of two patches of cells. Patch 1 consists of cells 1 through j and 
Patch 2 contains cells j + 1 and beyond where of course, j is a positive integer less than m, the total 
number of cells. The cells in Patch 1 may have different APD and/or CT restitution parameters 
than those in Patch 2. For simulations, we also assume that the pacer has a variable-rate within 
the aforementioned nominal sinus range. Assume that the pacer's oscillation period varies in a 
sinusoidal fashion: 

B(k) = Bq + Bx cos 

This beat pattern may represent an idealization (without random effects) of a single cycle of a 
complex, multi-pattern stretch of pacer beats. More complex variations may be considered in later 
studies if needed. 

The number w is the period of a full cycle of variable-rate oscillations, which we refer to as a 
full pacing cycle (FPC). The halfway point of the FPC, or the "bottom of the FPC well" occurs 
at k = w/2; at this point, the pacer oscillates with minimum period (fastest beat rate). On the 
other hand, when k = or k = w the pacer has its largest oscillation period (slowest beat rate). 
The number B\ is the FPC "amplitude" since it modulates the pacer's oscillation period within 
the FPC. If B\ = then the period is a fixed -Bo as in the preceding case-studies. 

Run 4.1. We fix the APD and CT parameters for Patch 2 to be the same as those used in 
Case-study 1. Patch 1 APD and CT have the same parameters as Patch 2 except for the following: 
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Cii 


0.2 


c 


300 ms 


170 ms 


0.05 cm/ms 



Thus APD is reduced (by a shift in the restitution function) in Patch 1; also the cells in this 
patch conduct more slowly with a lower value of c (the maximum CV). Additional loop and pacer 
parameters are set as follows: 



L 


B 


Bi 


w 


3 


5 


DIi 7 o, 1 < i < 135 


13.5 cm 


800 ms 


350 ms 


1000 


50 


140 ms 


600 ms 



The integer j gives the number of cells in Patch 1 and the three pacing parameters w, Bq and 
B\ define a variable pacing protocol with FPC period of 1000 beats and amplitude 350 about a 
fixed, nominal sinus trend of 800 ms. 

(i) The mode sequence of the SITR pattern obtained is 

{-464, [51], -1, [51], -1, [51], -1, [51], -1, [51], -2, oo}. 

(ii) The mode sequence exhibits a pattern of 5 incomplete reentry initiations followed by a 
pacer beat (or 2 in one case). Each incomplete reentry start which is blocked at the patch junction 
(Cell 51) corresponds to a short cycle length (about 224 ms) compared to the pacer's period of 
approximately 460 ms near the middle of the FPC. Hence the pattern of 5 alternating reentry/pacer 
beats resembles a series of premature ventricular contractions (PVC) in a bigeminal form before 
reentry finally takes hold for good. 

(iii) We note that if j = 40 or fewer, then no reentry is initiated in this run. On the other hand, 
larger j than 50 sustain reentry more quickly by increasing the percentage of slow cells in the loop. 

Run 4.2. In the preceding run the pacer did not interfere with the reentrant circulation (TP 
did not fail) because 5 = 140 ms was large enough to block the pacer. We now show that decreasing 
5 to promote pacer interactions may cause the eventual termination of reentry in the preceding run. 
We make only one change in this run relative to Run 4.1 by setting 5 = 70 ms. 

(i) The mode sequence of the resulting SITR pattern is as follows: 

{-464, [51], -1, . . . , [51], -1, [51], -2, 1, -1, 11, -3, 1, -1,4, -1, 6, -3, 1, -1, 

N v ' 

4 times 

4, -1, 6, -3, 1, -1, 11, -1, 6, -3, [51], -1, . . . , [51], -1, [51], -374, . . .} 

V v ' 

38 times 

After several PVC's and reentrant bursts of different durations, the pattern locks into the 
paced mode until reentry is initiated in the next FPC. This mode sequence is sensitive to computer 
round-off errors; however, the varied nature of mode changes persists. Figure 9 gives a graphical 
representation of mode changes for this run. 

(ii) The first deviation from the mode sequence of Run 4.1 occurs with a TP failure that 
interrupts reentry after just one full cycle. Some of the —1 entries in the mode sequence above 
(mostly the ones not associated with conduction blocks in Cell 51) indicate failures of TP also. 
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Figure 9: Schematic diagram of mode changes in Run 4.2. Beats in the reentry mode are listed with value 
1 and those in the paced mode are given the value 0. 

4 Further details of the model. 

In this section we present additional supporting details about the model. 

4.1 Remarks concerning the APD restitution function. 

4.1.1 Memory and latency. 

Experimental measurements and other studies reveal the difficulty of fitting a single curve to the 
data, particularly at the lower end of the DI scale. To address this and related issues, refinements 
have been considered in the forms of memory and latency. Beat to beat memory has been studied 
in [13, 19, 21, 23, 28, 46, 56] and has been modeled in a number of ways. In many cases, modeling 
memory can be accomplished by using DI values from several earlier stimulations. Latency (roughly 
a brief period between the arrival of a pulse and the firing of the cell) has also been expressed as a 
restitution function of DI [12, 13]. If this effect is added to APD we obtain a composite restitution 
function that may no longer be monotonic for small DI values. For simplicity, we ignore both 
memory and latency in this paper; however, they can (and need to) be included in the model for 
greater refinement in future studies. 

4.1.2 Standard vs. dynamic protocols. 

The APD restitution function may be measured experimentally using two different pacing protocols: 
Dynamic or standard (S1-S2) [29, 32]. These protocols generally produce APD curves that are 
shifted or have different slopes. For a discussion of this issue regarding the APD restitution that 
we use in this paper see [33]. The fit in ([1]) is to the dynamic restitution data in [33]. 
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4.1.3 Remodeling. 

The APD values (and thus the restitution function) may be affected by the steady fast pacing of the 
ventricles [34] or of the atria [39] . The implied shift in the APD restitution can have a significant 
effect on the SITR patterns in certain borderline cases. For simplicity we do not consider such effects 
in this paper; however, a suitable modification of the APD restitution function that incorporates 
explicit time dependence can accommodate dynamically-induced shifts. 

4.1.4 APD non-monotonicity and spatio-temporal chaos. 

Experimental evidence is offered in [53] that the APD restitution curve may be non-monotonic 
for small DI values. Working with a APD restitution curve having a single local minimum at 
low DI, the same study demonstrates the occurrence of spatial bifurcations along the length of a 
periodically paced fiber of cardiac tissue. Although not involving a loop structure, the work in 
[53] uses restitution functions in a discrete model. A localized form of non-monotonicity (like a 
bump or dent) in APD may occur at any DI value and has interesting consequences for additional 
complexity in reentry mode. In [48] where a continuous model of the loop as a cable is used, this 
type of non-monotonicity is shown to be responsible for the appearance of spatio-temporal chaos in 
APD and other key variables over the length of the loop. The occurrence of spatio-temporal chaos 
on the loop may raise the level of unpredictability in SITR if the amplitudes of oscillations in DI 
are large enough to let the wavefront reach its tail (thus causing a threshold failure). 

4.2 Restitutions of ionic currents and conduction block. 

Various anatomical and functional mechanisms for UB in a loop are discussed in the model study 
[49]. These mechanisms include a properly timed premature stimulus delivered to a suitable location 
in the loop, inhomogeneities in the degree of cellular uncoupling and gap junction resistance, in 
membrane excitability and in fiber cross sectional area. The discussion in [49] points to time, space 
and voltage "windows of vulnerability" for the occurrence of unidirectional block and initiation of 
reentry. Of these, the time window is seen to be the easiest to measure and its value is related to 
the space and voltage windows via standard mathematical formulas. 

Each cell requires a certain amount of depolarizing current to fire an action potential. This 
activation threshold can be represented by a restitution function [27] or equivalently, the strength- 
interval curve [13] shifted to the left so that the origin represents the end of the ERP. In Figure 10 
this curve is shown as a decreasing restitution function I a . 

In principle, at any time following the end of the ERP a sufficiently large current will be able to 
elicit an action potential. However, in the absence of external electrical shocks, there is a limited 
amount of current I g available to depolarize the next cell. If a restitution of I g is given then the 
intersection point of I a and I g is the DI threshold or cut-off value DI* in the sense that if DI > DI* 
then action potential is fired and propagation occurs. The time interval DI* is similar in nature to 
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Figure 10: Activation and generated currents restitution curves illustrating the definition of Dl* and Dl** . 
In this figure, I a = if = 1^ (see text). 

the chronaxie [47] although the definition of Dl* requires another restitution curve in addition to 
I a . A monotonically increasing restitution function for I g might be appropriate; see Figure 10. 

For our purposes here we do not need explicit formulas for either I a or I g . The value Dl* may 
be chosen arbitrarily or estimated experimentally without reference to its potential sources. Values 
as high as 160 ms were used in [26]. In human cases the value of Dl* in vivo is likely to be in a 
much smaller range. We use the conservatively small value of about 15 ms in this study. 

An analog of the time window in [49] can be defined here using the restitution functions I a and 
I g . We define two sets of currents restitution functions for Cell 1 corresponding to the two possible 
directions in the loop: I„ , lj? and Ig, lj* for the "left" and the "right" directions. The asymmetry 
creates a UB window in Cell 1 as follows: Let > 1% and Ig < Ig with at least one of these 
inequalities strict. Let Dl* be the (unique) intersection point of l£ and Ig and let Dl** to be 
the intersection point of 1^ and Ig. Then Dl* < Dl** and the interval between these two values 
is where the UB occurs. For if DI\ n represents the Dl of Cell 1 in a cycle or beat n, then the 
conditions 

Dl* < DI hn < Dl** 

imply that propagation is possible in one direction because the threshold marked by Dl* is crossed 
but it is inhibited in the other direction where the Dl** threshold is not crossed. 

We emphasize that the unidirectional block as defined here is functional and may not occur 
even if the UB window exists (for anatomical or functional reasons). If the number DI\ n does not 
enter the the UB window (i.e. the interval from Dl* to Dl**) in a particular beat n then the UB 
does not materialize. This option adds another mode to the loop-pacer system, namely the BDC 
mode mentioned previously. In addition to TC and TP, the following threshold is associated wtih 
the BDC mode: 
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TUB: D/i in _i < DP* with beat n — 1 in the paced mode. 

This inequality defines the UB threshold to ensure that Cell 1 conducts unidirectionally. If beat 
n — 1 is in reentry mode then conduction in beat n is automatically unidirectional whether TUB 
holds or not, because the pulse is transmitted from a still active Cell m. To limit the variety of 
cases to consider, we assume in this paper that DP* is so large that the unidirectional block is 
always present in Cell 1 (permanent UB) so TUB always holds. 

5 Summary and discussion. 

In the preceding sections, we developed an interactive loop-pacer threshold model that is a part of 
a familiar cardiac anomaly. Using experimental data on restitution parameters from the existing 
literature for action potential duration and conduction time, we obtained various conditions and 
parameter values that would cause irregular patterns of spontaneous initiations and terminations 
of reentry. 

A major cause of unpredictable behavior is threshold bistability where one stable regime causes 
a failure of the circulation threshold TC and another does not. In this situation the initiation and 
termination of reentry depends on the location of the DI state vector in the m-dimensional state 
space, where the nonlinearity in propagation equations makes it difficult to track the state vector. 
Thus, it is quite difficult to forecast the long-term behavior of an evolving SITR pattern for a 
T-bistable loop. Another potential source of unpredictability is the interference by the pacer. If 
reentry occurs, then the pacer and the reentrant circuit are thrown out of phase so a phase resetting 
pacer pulse may unexpectedly end the reentrant circulation. 

A number of prior studies of a different nature from ours also find complex SITR patterns in 
simple settings, such as in sheets [3, 25], loops [44] and aggregates [35] of animal cardiac tissue. 
In [3] and [25] complex patterns in monolayers of growing cell cultures are linked to changing 
densities and to local inhomogeneities, respectively. In [35] bursting beat patterns are studied 
in cell aggregates under external pacing. Closer to the subject matter of this paper, in [44] a 
ring-shaped lab preparation is considered with two localized pacemakers. Threshold relations for 
handling the pacers in this loop appear in the mathematical discussion. In a different class of 
studies, namely, the modulated parasystole, a permanent ectopic rhythm coexists with the sinus 
or another competing rhythm. In such relatively simple systems, complex beat patterns arise that 
are highly sensitive to parameter changes [23, 28, 43]. Like the TP threshold in our model, a 
mathematical relation with a discontinuity is responsible for the complexity of temporal patterns 
in parasystole; see [23, 43]. 

We omitted many features from our discussion here that one might include in an enhanced 
version of this model. For example, UB may be treated functionally by refining the UB time window 
and considering bidirectional propagation in the loop, as noted above. Other enhancements to this 
model might include adding memory and latency to the APD restitution function. The inclusion 
of memory is especially important in clarifying the dynamic significance of delayed responses (not 
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to mention improved fits to the APD restitution data). Additional directions for exploring the 
dynamics of the loop-pacer system include consideration of multiple exit / entry points into the loop 
that can affect reentry, a systematic approach to heterogeneity and also some theory to better 
model the time lag S n that influences the loop-pacer interactions. A somewhat different direction 
that may be explored would be to incorporate a programmed shock routine to artificially end long 
bursts of fast, reentrant beats. 

With some of the above enhancements the SITR patterns will likely become more complex and 
less predictable, though it is also possible that some of these enhancements (e.g. memory) can 
have a moderating effect. Understanding the precise mechanisms that generate such complex SITR 
patterns may in turn offer new insights into the nature and causes of long-term temporal patterns 
of tachyarrhythmia occurrences. 
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